En este notebook empezaremos con un micro-curso de probabilidad usando nuestra nueva herramienta. En este momento es importante recalcar que a medida que vayamos avanzando en los notebooks, las contribuciones en código de los autores serán cada vez menores, puesto que consideramos que lo aprendido en la primera parte del curso, más la sencillez de PyCUDA proveen al lector las herramientas necesarias para poder hacer sus propios códigos desde 0.
Al lanzar una moneda, uno puede decir con absoluta certeza que sólo una de las dos caras de la moneda verá hacia el frente, e incluso ir más allá: la probabilidad de caiga en cualquiera de las dos caras es de $\frac{1}{2}$.
Recomendamos leer sobre la visión "frecuentista" y "bayesiana" de las probabilidades que se le puede adjudicar a un evento. (Revisar las referencias)
Definición: Llamamos Distribución de Probabilidad a la colección de todas las probabilidades de todos los eventos de una variable aleatorio $X$ asociada a un evento.
Ejemplo: Al lanzar una moneda obtenemos dos posibles eventos.
De este manera obtenemos que la Distribución de Probabilidad asociada al lanzar la moneda es: $$\{ P(X = Águila) = \frac{1}{2}, \ P(X=Sol)= \frac{1}{2} \}$$
En este ejemplo, las probabilidades son iguales en cada valor de $X$, por lo que llamamos a esta distribución uniforme. Pero, ¿qué tan uniforme es?
Cambiaremos ahora a otro ejemplo un poco más complicado: el lanzamiento de un dado de 6 caras. No es difícil llegar a que, si el dado no está trucado, cada una de las caras tiene una probabilidad de $\frac{1}{6}$ de quedar mirando hacia arriba.
Pasemos ahora a la programación.
Como primer objetivo haremos que $N$ threads generen un número entero aleatorio entre $1$ y $6$.
In [25]:
import pycuda.autoinit
import pycuda.driver as cuda
import numpy as np
from pycuda.curandom import rand as curand
from pycuda.compiler import SourceModule
import pycuda.gpuarray as gpuarray
# Escribimos el kernel que utilizaremos
# El kernel toma un arreglo con N números aleatorios entre 0 y 1 y los convierte en números aleatorios
# entre 1 y 6.
mod = SourceModule("""
#include<stdio.h>
#include<math.h>
__global__ void Dado( float * d_arregloIn, float * d_arregloOut) {
int idx = blockDim.x * blockIdx.x + threadIdx.x ;
d_arregloOut[idx] = ceil( 6*d_arregloIn[idx] ) ;
}
""")
# definimos el número de dados a lanzar
N = 50
# Iniciamos los arreglos en el GPU
# CuRAND es una librería especial para lanzar números aleatorios y hacer cálculos estadísticos
a_gpu = curand(N,) # Ver la nota al final del Notebook
b_gpu = gpuarray.zeros_like(a_gpu)
# Rescatamos la función para luego ejecutarla
Dado = mod.get_function("Dado")
# El kernel es lanzado de forma muy parecida a CUDA C
Dado(a_gpu, b_gpu, block = (N, 1, 1), grid = (1, 1, 1))
print b_gpu
# Finalmente liberamos la memoria
a_gpu.free()
b_gpu.free()
En este primer código creamos únicamente un lanzador de $N$ dados de $6$ caras. Ahora pasemos a calculos un poco más interesantes.
[1] Como primer ejercicio tendrás que modificar el código de arriba de tal forma que en el GPU se calculen las frecuencias de cada uno de los valores obtenidos. El resultado final tiene que ser un arreglo en el cual estén las frecuencias.
Puedes basarte en el siguiente código de Python
F = np.zeros(6)
N = 50
for i in resultado:
for j in xrange(1,7):
if i == j:
F[i] += 1./N
Puedes graficar estos valores en un histograma usando la funcion
np.hist()
para luego repetir el experimento con el mismo valor de $N$ y luego aumentando el valor.
¿Qué observas?
[2] En el histograma del ejercicio pasado se logra observar cómo entre mayor es el número de tiros, más se acerca la frecuencia de cada cara a la frecuencia teórica ($\frac{1}{6}$). Ahora tu trabajo será de calcular la desviación entre la frecuencia obtenida y el valor teórico. Evidentemente este cálculo debe de hacerse en el GPU.
Ahora buscamos observar como graficar observar la desviación entre la frecuencia obtenida y su valor teórico. Al graficar directamente la desviación obtenemos una gráfica en la que no se observa mucho puesto que hay un decaimiento muy rápido para luego converger hacia 0. Sin embargo, al realizar una gráfica log-log obtenemos una recta tal como se muestra en la figura siguiente. ¿Qué significa esto? La respuesta es que la desviación estándar se comporta como una potencia $\sigma(n)=n^{-p}$.
Esta potencia es $p=-\frac{1}{2}$ y es debido a que la simulación que realizamos es Monte Carlo. Todas aquellas simulaciones de esta tipo tendrán una desviación estándar que tendrá este tipo de comportamiento.
(a) Pensemos sólo en el lado $1$. Sea $p := \frac{1}{6}$ la probabilidad de que caiga en este lado. Si lanzamos 2 veces, ¿cuál es la probabilidad de que salgan (i) cero veces $1$; (ii) una vez $1$; (iii) dos veces $1$?
(b) Si lanzamos 3 veces, ¿cuáles son las probabilidades de que salga cada número posible de $1$s entre 0 y 3?
(c) Si lanzamos N veces, ¿cuáles son las probabilidades?
(d) ¿Cuál es el valor promedio que sale? ¿Qué otra cantidad podríamos querer calcular?
Ahora pasemos a otro tipo de distribución de probabilidades basada en el ensayo de Bernoulli. La idea principal es sencilla: Un cierto evento tiene una probabilidad $p$ de éxito. Si el evento es exitoso, entonces se marca con el valor $1$. Sino, se obtiene el valor $0$.
Tomemos el ejemplo del dado. El evento exitoso será el obtener la cara con valor $1$, lo cual quiere decir que tenemos una probabilidad $p = \frac{1}{6}$ de que nuestro evento sea exitoso.
Al repetir el lanzamiento $N$ veces, entonces podemos obtener la distribución de probabilidad $Y_{N}$, la cual nos dará la colección de todas las probabilidades de que el evento sea exitoso cuando se repite $i$ veces el experimento. $i = 1, ..., N$
Para todo lo siguiente supondremos que nuestro evento exitoso será obtener un $6$ a la hora de tirar un dado.
[3] Ahora pasemos a algo un poco más complicado. Escribe un kernel que haga $N$ ensayos de Bernoulli y que calcule la proporción de eventos exitosos.
Para este problema tendremos distintas manera de obetener el resultado dentro del GPU(en el CPU el problema es trivial).
La manera fácil y no eficiente es escribir un kernel en el cual sumemos un $1$ si el evento fue exitoso. Esto puede hacerse con un condicional if.
En el caso en el cual nos decantemos por algo más eficiente, tendremos que recurrir a un kernel más sofisticado. El tipo de operación que se hará se llama "reducción" y ha sido ampliamente discutido a lo largo de blogs y foros.
La idea es esta: supongamos en un primer tiempo que tenemos un bloque de 8 entradas. Entonces el procedimiento será sumar por pares -tal como se muestra en la fiugra siguiente- hasta poder obtener el resultado final del bloque.
Ahora, si tenemos un arreglo de 8 bloques entonces el siguiente paso será llevar cada uno de los resultados de los 8 bloques a uno solo. De esta manera obtendremos un bloque con los resultados de los demás bloques.
El paso siguiente será de nuevo, realizar el procedimiento de la figura anterior para obtener el resultado final.
Si el lector desea implementar y aprender esta herramienta, en las referencias colocamos un documento de NVIDIA en el cual se explica paso a paso un kernel de reducción muy sofisticado. Esta alternativa es altamente recomendable por todos los usos que se le puede dar.
[4] Finalmente calcula la distribución de probabilidad $Y_N$ con los resultados obtenidos en el inciso anterior.
Aquí te recomendamos usar distintos kernels, uno __device__ con el que calcules la proporción de eventos exitosos, y otro __global__ con el que calcules $Y_N$.
La distribución de probabilidad de $Y_{N,p}$ se llama una distribución binomial. Se dice que es una distribución de probabilidad discreta, ya que los valores que toma son valores discretos.
La notación que se utiliza es
$$X \sim B(N,p)$$para decir que $X$ es una variable aleatoria con distribución de probabilidad binomial con parámetros $N$ y $p$.
En este primer Notebook, y a partir de ahora usaremos la librería cuRAND la cual nos permite generar números aleatorios para cada thread dentro del GPU. Esto tiene grandes ventajas con respecto a la utilización de la función rand() dentro de numpy, sobre todo debido al mayor control que nos da cuRAND sobre los generadores. Sin embargo, cuRAND tiene la desventaja que el máximo número de generadores activos en GPU's tipo Tesla es de 256. Para el caso de GPU's tipo Fermi, 1024 generadores pueden ser activados sin -en principio- tener mayor problema.
Recomendamos ampliamente leer las dos primeras partes del Apéndice 3 para aprender un poco más sobre como se generan números pseudo-aleatorios y cuasi-aleatorios en el GPU con ayuda de cuRAND.
In [ ]: